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ABSTRACT 

We introduce a sub-grid force correction term to better model the dynamical fric¬ 
tion (DF) experienced by a supermassive black hole (SMBH) as it orbits within its host 
galaxy. This new approach accurately follows a SMBH’s orbital decay and drastically 
improves over commonly used ‘advection’ methods. The force correction introduced 
here naturally scales with the force resolution of the simulation and converges as 
resolution is increased. In controlled experiments we show how the orbital decay of 
the SMBH closely follows analytical predictions when particle masses are significantly 
smaller than that of the SMBH. In a cosmological simulation of the assembly of a small 
galaxy, we show how our method allows for realistic black hole orbits. This approach 
overcomes the limitations of the advection scheme, where black holes are rapidly and 
artificially pushed toward the halo center and then forced to merge, regardless of their 
orbits. We find that SMBHs from merging dwarf galaxies can spend significant time 
away from the center of the remnant galaxy. Improving the modeling of SMBH orbital 
decay will help in making robust predictions of the growth, detectability, and merger 
rates of SMBHs, especially at low galaxy masses or at high redshift. 

Key words: Numerical methods: Supermassive black holes: cosmological simulations: 
dynamics 


1 INTRODUCTION 

Supermassive Black Holes (SMBHs) represent a crucial, 
though still poorly understood, aspect of galaxy evolu¬ 
tion theory. SMBHs with as much as 1 0^ Mq in mass 
power luminous quas ars observed at z > 6 llFan et al.ll200ll : 
iMortlock et al.lboill l. while in the local Uni verse, SMBHs 


are UDiquitous in Dotn massive gaiaxies i e.e. iGrenren 
1984 iKormendv & Richstondll995l: iKormendv & Hoi 

et ai.i 

2013) 

and 

small, bulge-less disk galaxies (Shields 

et al. 

2008; 

Filinoenko & Hoi l2003h as well as dwarfs ( 

Reines et al. 

2011 

: Reines & Dellei 20121: iReines et al.ll2013l: 

Moran et al. 

2014 

SME 

host 

). Empirical scaling relations between the mass of 
)Hs and the stellar mass and velocity dispersion of their 
galaxies are indicative of co-evolution (iHaring & Rix 

20041; Giiltekin & et al.ll2009l: ISchramm & SilvermanI 20131: 


Kormendv &: Hoi l2013l l. Powerful outflows resulting from 
feedback from an active nucleus have been observed 
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llFeruglio et al.l 120101 : lAlatalo et al.l l201lh and are thought 
to play a critical role in the quenching of star formation in 
massive galaxies and the exponentia l cutoff of the galaxy 
luminosity function at high masses llSpringel et al. l2005al : 
iBower et al.ll200'^ : ICrotonlT2009l : iTevssier et al.ll2011 ). 


Numerical simulations have proven to be a critical tool 
for understanding the formation and evolution of galaxies 
and extensive work has already be en done to implement 
SMBHs into these simulations (e.g. Di Matteo et aDl2003l : 
iHopkins et al.l[2005l : ISiiacki et al.ll2009ll . Due to the scale of 
these simulations, which often include a full cosmological 
environment, the resolution is necessarily limited. Accretion 
onto SMBHs and the following feedback processes are hence 
implemented via sub-grid prescriptions, under the broad as¬ 
sumption that conditions at the smallest resolved scale drive 
the BH evolution at the scale of just a few parsecs. One ma¬ 
jor obstacle is that it becomes necessary to accurately follow 
the dynamics of a single object (the SMBH), something that 
cosmological simulations with force resolution larger than a 
few pc are inherently not well equipped to do. These dy- 
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namics can have important consequences for the growth of 
SMBHs, particularly during galaxy mergers. Detailed simu¬ 
lations of galaxy mergers show that the dynamics of SMBHs 
are intimately conne cted with accretion and vice versa (e.g. 
ICallegari et al.lboilf l. 

In this paper, we describe our approach to take into ac¬ 
count unresolved dynamical friction (DF) in the orbital evo¬ 
lution of SMBHs and show that our prescription is a promis¬ 
ing and realistic alternative to advection that can easily be 
implemented in existing codes. We show that this approach, 
when applied to simulations with the typical resolution of 
modern cosmological simulations 0.1 — 0.5 kpc) leads to 
more realistic dynamics that match well with analytic ap¬ 
proximations for DF timescales. In Section 2 we summarize 
some current methods for correcting SMBH dynamics in cos¬ 
mological simulations and why a new method is needed. In 
section 3 we present our method. In section 4 we present 
our results, comparing our method with a popular advec¬ 
tion prescription in both an isolated dark matter (DM) halo 
and a fully cosmological zoomed-in simulation of a dwarf 
galaxy. In section 5 we summarize and discuss our results. 


2 THE NEED FOR A NEW MODEL 


Dynamical friction, the force exerted by the gravitational 
wake cau sed by a massive obj e ct moving in an extended 
medium (IChandrasekhail 119431 : iBinnev fc Tremaind l2008l l 
causes the orbits of SMBHs to decay towards the center of 
mass ive galaxies llGovernato et all 1 19941 : iKazantzidis et ahl 
I 2 OO 5 II . In cosmological simulations, modeling DF is partic- 
ulary challenging, as the particle representing a SMBH is 
often only a few times the mass of the background par¬ 
ticles. In this case SMBH orbits can become dynamically 
heated due to the limited mass resolution and the result¬ 
ing spurious collisionality from a noisy gravitational poten¬ 
tial llHernauist fc Barnedll99^ . This numerical heating can 
cause the black holes to gain/lose energy and be unrealis¬ 
tically perturbed away from the center of its host halo. In 
order to lessen collisionality and improve performance, N- 
body simulations employ a gravitational softening length, 
a characteristic scale below which the gravitational force 
between two particles becomes damped. This mechanism, 
while necessary numerically, hinders DF by preventing the 
close interactions that are necessary to form a wake in the 
vicinity of the SMBH. 

To account for the unrealistic dynamics of SMBHs 
in cosmological simulations, many groups employ artih- 
cial advection schemes that change the motion of SMBHs. 
Each method comes with its own drawbacks. For exam¬ 
ple, placing the SMBH at the posit ion of the lowest po¬ 
tential gas part i cle a round it (e.g. ISpringel et al.l l2005bl : 
iBooth &; Schavd l2009l ~l causes chaotic motions, especially 
when relative velocity constrain ts are put on the gas par¬ 
ticles (IWurster fc Thackeil I20l3 l. Utilizing a large ‘tracer 
mass’ with which the SMBH gravitat ionally interacts with 
its surroundings llPebuhr et al.ll20Ill l affects the morphol¬ 
ogy in the g alactic center as well as th e accretion history of 
the SMBH dWurster fc Thackeil l2013l ~l. Pushing the SMBH 
a certain distance in the direction of ei ther the local stellar 
density g radient dOkamoto et al.l l2008tl or the local center 
of mass dWurster &: Thackeil 2013ll can avoid chaotic mo¬ 


tions but adds an additional free parameter, the distance the 
SMBH is pushed each step. Even disregarding their individ¬ 
ual drawbacks, none of these solutions are ideal, as they all 
require the explicit assumption that SMBHs rapidly decay 
into and remain stable at th e center of their host galaxies , 
which is often not accurate dBellovarv et al.ll20ld) . 

In strongly interacting systems, where SMBHs are 
thought to go throug h much of their g rowth and pos¬ 
sibly their formation dMaver et al.l 120071 ). the inner re¬ 
gions of a galaxy are likely to experience strong poten¬ 
tial fluctuations, affecting the dynamics of the SMBH. This 
could be es pecially evident in dwarf galaxies, whe re black 
holes exist dReines et al.l l20I3l : iMoran et al.l I20l4l within 
a shallow potential with a cored density profile I Oh et ahl 
I 2 OIII : iTevssier et al.l[20I3h . making perturbations to their 
orbits more likely. Additionally, SMBHs accreted during 
both major and minor mergers, especially dry mergers 
dKazantzidis et al]l2005l l. do not immediately sink to the 
center of the ga laxy, creating a population of small ‘wan¬ 
dering’ SMBHs dlslam et al.112003 : IVolonteri fc Pern3l2005l : 
iBellovarv et al.ll2010h . Advection, by assuming that orbital 
decay timescales are always short, could then artificially in¬ 
crease the mass of central SMBHs, through inflated merger 
and accretion rates. 

M ore advanced approa ches have been explored. For ex¬ 
ample, iDubois et al.l d2013l l include a sub-grid prescription 
for gas dynamical friction acting on the SMBH. While this is 
a promising method, it requires assumptions about the mul¬ 
tiphase nature and equation of state of the gas far below the 
resolution l i mit o f any cosmological simulation (~ 1 — 5 pc). 
iLupi et ni d2015fl increase the resolution around SMBHs to 
attain very accurate dynamics, but their method is applica¬ 
ble only for very high resolution simulations meant to closely 
follow BH-BH mergers, not cosmological simulations. 

The method we propose in this Letter is a simple solu¬ 
tion to correcting SMBH dynamics in cosmological simula¬ 
tions. Our approach is to estimate the unresolved dynamical 
friction felt by the SMBH and apply the appropriate force to 
the SMBH particle. This is a significant improvement over 
artificial advection, as it makes no assumption about where 
the SMBHs should be located in a galaxy at a given time 
and it has no explicit effect on the SMBH surroundings. In 
addition, it requires minimal assumptions about the state 
of the simulation below the resolution limit and it naturally 
converges with increasing resolution. 


3 METHODOLOGY 
3.1 The Test Simulations 

To test the dynamics of black holes in a realistic DM halo, 
we generate initial conditi ons for a collaps ing overdensity 
following the procedure of lEvrar 3 lll988l l^. We begin the 
simulation by approximating the state of a halo beginning to 
collapse at high redshift. We then allow the halo to collapse 
before implanting a black hole of mass 10® Mq. We run a 
suite of simulations at different mass and spatial resolutions 


^ Created using the package ICInG created by M. Tremmel. 
The code is publicly available and can be downloaded at 
https: //github.com/mtremmel/ICInG.git 
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Table 1. Information about the resolution of each isolated dark 
matter halo test simulation. 


Name 

Number 

of Particles 

Dark Matter 
Mass (Mq 

Softening 
length, Eg (pc)"^ 

Low Res 

1.05 X 10® 

9.78 X 10® 

311 

Oversampled 

8.39 X 10® 

1.22 X 10® 

311 

High Res 

6.71 X 10 '^ 

1.53 X 10"^ 

77 

Small Soft 

6.71 X 10’’ 

1.53 X 10"^ 

10 


* Mass of the dark matter particles in the simulation 
t spline kernel gravitational softening 


(see Table 1) with a massive ‘black hole’ particle initially 1) 
at the center while the halo is still actively collapsing or 2) 
on an eccentric orbit after the halo has mostly relaxed. The 
resolution of the ‘Low Res’ and ‘Oversampled’ r uns are com¬ 
parable to current uniform volume simulations llKeres et al 


20]^ [^iiacki et akl l2014l : ISchave et al.l l2015l : iDubois et al 


2014al L while the ‘High Res’ simulation is representative 


of high resolution ‘zoomed-in’ s i mulations llGove rnato e t al.l 
|2012 k IChristensen et al] l2014l : iHopkins et al. 20141 ). The 
force softening lengths adopted are typical of cosmological 
runs, while in ‘Small Soft’, a variant of the ‘High Res’ simula¬ 
tions, the force resolution i s only lOpc, typical of simulations 
of isolated binary mergers llCanelo et al.ll2015l) or very high- 
res co smological zoomed-in simulations (e.g. iDubois et al.l 
l2014bl) . 


We verified that the halo formed in the collapse has a 
density profile typical of CDM halos llNavarro et al.|[i996l) . 
The halo we use here has (after a time consistent with z ~ 0) 
yi^ir ~ 2 X 10^^ Mq, R„i,. ~ 115 kpc, and a concentration 
c ~ 4.5. An analytic expression of the approximate timescale 
for a rigid object to sink to the center of a DM halo was cal¬ 
culated bv lTaffoni et al.l ll2003l) . Given these conditions and 
the mass of the test black hole, the estimated DF timescale 
(tdf) for an eccentric (v= O.lvcirc) orbit at an initial dis¬ 
tance of 2 kpc from the halo center is approximately 1.8 
Gyr. This initial orbit was chosen as typical of a SMBH af¬ 
ter its parent satellite has been tidally disrupted. It has a 
non-negligible tdf, but still much less than a Hubble time. 


3.2 The Dynamical Friction Prescription 

Dyna mical friction occurs due to both large (iGolni et al.l 
I 1999 II and small scale perturbations to the black hole’s sur¬ 
roundings. We consider perturbations on scales larger than 
the gravitational softening length, eg, to be well resolved. 
On these scales the potential should be smooth, so long as 
enough particles are used so that dynamical heating is min¬ 
imized. 

The acceleration a black hole of mass Mbh feels from 
DF due to particles of mass ma <S Mbh is given by the 
following formula: 


SLDF = -47rG^MsirmolnA / T■Vaf{'Va)-r—^ -(1) 

J IVBff— VaF 

Where /() is the velocity distribution function, InA is 
the Coulomb logarithm, Va is the velocity of the surround¬ 
ing background objects relative to the local center of mass 
(COM) velocity, and vbb is the relative velocity of the black 


hole. We assume that within tg from the black hole the ve¬ 
locity distribution is isotr opic, giving Chandras ekhar’s Dy¬ 
namical Friction Formula llChandrasekhailll943h . 

SLDF = -47rG^MmalnA^^^ f dvavlf{^ra) (2) 
"^BH Jo 

This can be further simplified by substituting the inte¬ 
gral for p{< vbh), which is the density of particles moving 
slower than the black hole. 


bldf = —47rG^Mp(< VBB)lnA (3) 

^BH 

The Coulomb logarithm depends on the maximum and 
minimum impact parameters, hmax and hmin, such that 
InA ~ ). Because DF is well resolved at scales 

greater than the softening length, we set hmax = tg to avoid 
double counting frictional forces that are already occurring. 
For the minimum impact parameter, we take it to be the 
minimum 90° deflection radius, with a lower limit set to the 
Schwarzschild Radius, RscJi- 


bmin = maa:(b9o,Rsc/i);b9o = —5 - (4) 

VSB 

For the calculation, we use 64 collisionless particles (i.e. 
dark matter and star particles, if present) closest to the black 
hole. We calculate the velocity of each particle relative to the 
COM velocity of those 64 particles. We verified that our re¬ 
sults do not depend strongly on the number of neighbors 
used, although using too few particles could result in nu¬ 
merical noise in the calculation of this force. Since we are 
explicitly assuming the velocity distribution is isotropic, the 
following must be true. 


p(< vbb) 


M(< Vbb) 

Mtotal ^ 


(5) 


Where p, the total density around the black hole, is 
calculated by smoothing over the chosen 64 particles, i.e. 
p = rriiW (jbh — Vi, h). M(< vbb) is the total mass of 
the chosen particles that are moving slower than the black 
hole relative to the local center of mass velocity and Mtotal 
is the summed mass of all 64 particles. 

The resulting acceleration (from eq. 3) is added to the 
black hole’s current acceleration, to be integrated the follow¬ 
ing time step. As the spatial resolution or black hole mass 
increases (or the velocity of the black hole decreases) hmin 
will become greater than hmax, in which case we claim DF 
is being fully resolved and therefore the correction is not 
needed. 

This method is not accounting for DF from gas, which 
can have important effects for supersonic black holes in re- 
gions where gas density dominates stars and dark matter 
llOstrikeilll999l : Tchapon et al.ll2013l) . This should not occur 
often on the scales relevant in these simulations. The center 
of larger galaxies are dominated by stars and smaller galax¬ 
ies have significant star and dark matt er fractions within 
the central regions dOh et al.ll2008l . l201ll ). so this effect will 
only be a minor correction in most cases. DF may be over¬ 
estimated within resonant DM cores where DF can become 
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No Added Correction 


With DF Correction 



Time (Gyr) 


Figure 1. Effects of the DF correction: Distance of a test 10® Mq black hole from halo center as a function of time at different 
resolutions. Dashed-dot and dashed lines indicate tg and 2eg respectively for both the Low Res and Oversampled models. The black hole 
starts on an eccentric (v= 0.1 Vdrc) orbit with apocenter of 2 kpc. The vertical solid line represents the analytically derived timescale 
for orbital decay. When the DF correction is applied, marked improvement is seen for all models except ‘Low Res’, which experiences 
too much dynamical heating due to the lower mass resolution. The orbits of ‘High Res’ and ‘Small Soft’ are very nearly the same once 
the correction is implemented, indicating numerical convergence when the DM particle mass is ~ 10^ Mq. 


much less efficient (|R.ead et al.l 1200(^ 1 . This effect is sec¬ 
ondary and mainly important when the gravitational soft¬ 
ening length is appreciable compared to the size of the core 
structure. Often the orbital differences should be smaller 
than or similar to the resolution limit and therefore unim¬ 
portant. Additionally, interactions with clumpy gas has been 
shown to significantly increase the timesc ale for the or¬ 
bital decay of SMBH bin aries below ~ lOOpc llFiacconi et al.l 
I2OI3I : iRoskar et al.ll2015l l. However, this effect would not be 
well resolved by even high resolution ‘zoomed-in’ cosmolog¬ 
ical simulations. 

4 RESULTS 

4.1 Isolated Dark Matter Halo 

We find that, for SMBHs placed in the center of a collaps¬ 
ing halo, only the Low Res simulation experiences significant 
dynamical heating, causing the BH to be unrealistically per¬ 
turbed away far from halo center. In the rest of the runs the 
BH remains within one softening length from halo center 
for 6 Gyrs and shows no sign of heating. This is due to the 
higher mass resolution present in those runs. Simulations in¬ 
cluding either our DF correction or an advection correction 
(see below) have the same results, showing little difference 
from simulations with no correction at all for this scenario. 

Figure 1 shows the orbital evolution of a black hole ini¬ 
tially on an eccentric (v= 0.1 Vdrc) orbit with apocenter of 


2 kpc, placed in the halo after it has finished most of its col¬ 
lapse. The center of the hal o is defined at each step using the 
shrinking spheres method dPower et al.l 1200^1 . We verified 
that the density maximum and potential minimum coincide 
within much less than the force resolution. The vertical line 
represents the dynamical friction timescale for the orb it de¬ 
rived from the analytic model of iTaffoni et al.l ll2003tl and 
the horizontal lines represent eg and 2eg. Without the DF 
correction, only the ‘Small Soft’ model, with 10 pc spatial 
resolution and DM particle mass almost 100 times smaller 
than the BH, is able to show substantial orbital decay within 
2tdf- Implementing our DF correction results in a notice¬ 
able improvement for the orbital decay, even at the relatively 
modest resolution of the Oversampled model, where it falls 
to within 2eg of the center before 2tdf- At higher resolution, 
the dynamics converge to closely match with the analytical 
approximation. Note that even for our highest resolution 
simulation. Small Soft, the DF correction causes the SMBH 
to sink almost 1 Gyr sooner. These are very encouraging re¬ 
sults, as they Indicate that this correction results in realistic 
black hole orbital evolution even at resolutions attainable in 
large volume simulations and it has important consequences 
even at the highest resolutions tested here. 

Figure 2 compares the performance of our DF prescrip- 
tion to that of a co mmonly used advection scheme used in 
ISiiacki et al.l ll2007tl and various other simulations. The test 
is done using the Oversampled run, as this most closely re¬ 
sembles the resolution of a cosmological volume simulation. 
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Figure 2. DF correction vs Advection: Results from 
the ‘Oversampled/ model when implementing our DF correction 
(blue) compared with a commonly used advection routine (red) 
and no correction to dynamics (cyan). Dashed-dot and dashed 
lines indicate eg and 2eg respectively.The black hole starts on an 
eccentric (v= 0.1 Vdrc) orbit with apocenter of 2 kpc. Advection 
immediately pushes the off-center black hole to the center, miss¬ 
ing the orbital decay that our method captures well. Without any 
correction, the orbit decays far too slowly, remaining far (> 2eg) 
from halo center even after 6 Gyr. 

The BH is placed on an eccentric orbit, as in Figure 1. The 
advection scheme adopted repositions the black hole each 
time step to the position of the lowest potential particle 
within its 32 nearest neighbors while keeping the velocity 
unchanged. Not surprisingly, this results in the black hole 
staying very close to the center of the halo even when ini¬ 
tially set on an off-center orbit. The DF correction captures 
the more gradual orbital decay that the advection scheme 
completely misses. The run with no dynamical correction 
fails to have the BH sink within 2eg even after 6 Gyr. 

This is an important conclusion because these different 
orbital evolutions would result in drastically different ac¬ 
cretion histories for the black hole. Off center BHs should 
accrete less due to lower gas densities. Additionally, simu¬ 
lations that utilize advection would have black holes merge 
much sooner than what is predicted by their orbital decay 
timescale. Our improved method should then have impor¬ 
tant implications for the growth and merger rate of SMBHs 
in cosmological simulations of galaxy formation. 

4.2 Cosmological Dwarf Galaxy Simulation 

As a first test of the dynamics of SMBHs in a fully cosmo¬ 
logical setting, we run a high resolution ‘zoomed-in’ simula¬ 
tion that results in two dwarf galaxies with masses ~ 10^° 
Mq at z = 0. The simulation has a resolution similar to 
our High Res isolated halo model, with dark matter par¬ 
ticle mass 1.6 x 10"* Mq, gas particle mass 3.3 x 10^ Mq, 
and gravitational softening of only 87 pc. We showed in the 
previous section that at this resolution our DF prescription 
gives results that match analytic models. We chose a dwarf 


5 . 04.0 3.0 2.0 1.0 0.5 



5 . 04.0 3.0 2.0 1.0 0.5 



Figure 3. DF correction in a cosmological dwarf sim¬ 
ulation: The dynamics of four black holes in the cosmological 
zoomed-in dwarf galaxy simulation with DF (top) and advection 
(bottom). These are the black holes that end up in the most mas¬ 
sive system by the end of the simulation. Each colored line traces 
the distance of a black hole from the center of the most mas¬ 
sive halo. Black dots mark merger events and the dashed lines 
mark the gravitational softening length of the simulation (87 pc). 
Which of the two black holes emerges from a merger event and 
which is ‘eaten’ is unimportant. DF is able to sustain a long-lived 
dual black hole system (blue and red) while the advection scheme 
causes them to quickly merge. The green black hole remains on a 
very wide orbit in the DF run, but is quickly and unrealistically 
pulled to the center with advection. 


galaxy for this test because SMBHs are more likely to be¬ 
come perturbed away from galactic center, given their shal¬ 
low gra vitational potential and actively evolving, cored DM 
profile llGovernato et al1l2012l : IPontzen fc Governatdl20l3) . 
This will guarantee a useful test environment for exploring 
the differences between out method and advection. 

This test is also topical, as there is a growing sam¬ 
ple of dwarf galaxies w ith detected SMBHs llReines et al.l 
I 2 OI 3 I : iMoran et al.ll20i3l . Realistic numerical studies of BH 
formation and growth in these small galaxies, focusing on 
their occupation fraction and how they and their host galax¬ 
ies evolve toward the correlation with the stellar veloc¬ 
ity dispersion, would provide vital co nstraints on BH seed 
masses and early growth mechanisms IVolonteri et al.ll2008l : 
IVolonteri fc Gnedinll200^ : IVolonterilbOld) . 

We use the new N-body + SPH code ChaNGa 
dMenon et al.|[2015l ). which includes all of th e physics mod¬ 
ules previously Implemented In Gasoline dWadslev et al.l 
l2004h such as hydrodynamics, gas cooling, a cosmic UV 
background, star formation and SNe feedback. The ‘zoomed- 
in’ approach preserves the large scale tidal field while allow- 
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Figure 4. BHs in a Cosmological Dwarf: A snapshot of a zoomed-in cosmological simulation of a forming dwarf galaxy at z 
= 0.846. The gas density integrated along the line of sight is shown with darker colors indicating higher densities. In the dynamical 
friction simulation, a previous merger has created a central binary BH system (red and blue, see Figure 3). The separation of ~ 1 kpc 
is well resolved by the simulation, which has a force resolution of 87 pc. A more recent merger has set a third BH (in green) on a wide 
orbit (see Figure 3). In the image, the green BH is at its closest approach to galactic center. In the advection simulation, all BHs are 
quickly pushed to the center, where they merge, causing the simulation to miss these more realistic BH orbits. 



-10 -5 0 5 10-10 -5 0 5 10 


X[kpc] X[kpc] 


ing us to model a small region at hig h resolution. This issim- 
ilar to the simula t ion d escribed in iGovernato et al] ll2010l l 
and in IShen et alj ll2014l l. In this particular simulation, stars 
and DM densities dominate that of gas within the inner re¬ 
gions by more than a f actor of 10. Examinin g other dwarf 
galaxy simulations fe.g. IGovernato et alj|20foh . we find that 
gas densities can often reach similar values to DM and stars. 
In either case the contribution from gas to DF is only a mi¬ 
nor correction, at most a factor of a few, and negligible in 
our current example. 

By z <1 these systems are a good representation of real 
dwarf galaxies, with an extended stellar disk, no bulge, a 
high gas fraction and a cored DM profile. We first run the 
simulation until z = 6, when we insert five black holes of 
mass 5 x 10® Mq in the centers of the five most massive 
halos at the time, which all have a mass of 2 x 10® Mq 
or higher. In these simulations black holes do not accrete or 
produce feedback, as we are only interested in following their 
dynamics. From z=6 we run two simulations to z < 1, one 
with our DF routine and the other with advection. By the 
end of the simulation, there are only two major star forming 
galaxies with masses ~ lO*^® Mq. 

We then run the Amiga Halo Finder (AHF) 
llKnollmann &: Knebell200^ for all the saved snapshots and 
calculate the center of the main halo at each step using the 
shrinking spheres approach. In Figure 3 we follow the tra¬ 
jectories of four black holes with respect to the center of this 
halo (which originally just has the blue BH at its center). 
Each color represents a different black hole. Black dots in¬ 
dicate a black hole merger, which happens when two BHs 
come within two €g of one another at relative speeds low 


enough to be gravitationally bound. The dashed black lines 
indicate the gravitational softening length of 87 pc. 

The black holes become perturbed as the red, cyan, and 
blue host galaxies interact between z = 3 and 4. In the 
advection case the black holes are driven quickly toward the 
center where they all merge, leaving only one black hole 
(labeled as red). With the DF correction, only the cyan and 
blue BHs merge. After the red and blue galaxy hosts merge 
with DF, the blue and red BHs remain orbiting around the 
center of the merger remnant. The blue BH comes back to 
the center only after 4 Gyr and the red BH remains orbiting 
at around 1 kpc (11 times the force resolution) for another 
4 Gyr before sinking and merging with the blue BH. 

The more striking difference between the DF and the 
advection run involves the green black hole. When the much 
smaller green host galaxy merges with the blue/red host at 
2 ~ 1, it is initially far from halo center (~ 30 kpc) and 
is quickly disrupted by the main galaxy. With DF, the BH 
stays on a wide orbit, never coming much closer than a few 
kpc from halo center. In the advection case, however, the 
green black hole is quickly pushed to the center where it 
merges with the central red BH (see Figures 3 and 4). This 
is an unrealistic result, as the DF timescale of a 5 x 10® 
Mq black hole that far away from the center of such a small 
galaxy would be longer than a Hubble time. 

With this simulation we can clearly see how the choice 
of dynamical correction can affect the ability of SMBHs to 
become perturbed during mergers. In the DF simulation, 
BHs are able to remain off-center for many Gyr while with 
advection they are quickly driven to the center. Additionally, 
the DF simulation allows for sustained wide orbits resulting 
from minor mergers. Such dynamics can have an important 
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impact on interpreting the connection between the initial oc¬ 
cupation probability of SMBH seeds in dwarf gala:xy progen¬ 
itors and the observed occupation at low redshift. Methods 
such as the advection scheme presented here would predict a 
more direct connection, while the simulation with DF indi¬ 
cates that the nature of the mergers (i.e. mass ratio and ori¬ 
entation) can have an impact on which dwarf gala:xies have 
observable SMBH activity and when that activity occurs. 
The DF correction could also have important implications 
for BH merger rates, allowing them to become more decou¬ 
pled from galaxy merger rates than advection simulations. 
This will affect predictions of gravitational wave detections 
as well as estimates for recoiling BHs, which can have an 
impo rtant effect on observability (e.g. iMadau &: Quataerd 
|2004| . 

This simulation is a useful illustration of the variety of 
different, realistic BH orbits our method allows compared to 
commonly used advection schemes. In future work, we will 
explore the dynamics of BHs in a variety of different merger 
events within a cosmological context. 


5 DISCUSSION AND SUMMARY 

We have introduced a sub-grid force correction term for 
SMBH motion based on dynamical friction. This correction 
allows us to better model the orbital decay of SMBHs in 
numerical simulations. We have shown using controlled ex¬ 
periments of isolated DM halos that this addition matches 
analytic predictions of the orbital decay in DM halos with 
resolutions attainable by large-volume cosmological simu¬ 
lations. We have also demonstrated that our prescription 
naturally converges with resolution. 

This method is a significant improvement over exist¬ 
ing ‘advection’ methods that force a short orbital decay 
timescale regardless of the dynamical state of the system. 
When applied to a cosmological dwarf galaxy simulation, 
our method results in noticeably different black hole dy¬ 
namics compared with the advection scheme. In particular, 
our prescription: 

• Models the perturbation and gradual orbital decay 
of a central BH during a galaxy merger 

• Allows for long-lived dual BH systems with close 
(< 1 kpc) orbits. 

• Maintains a stable central BH when appropriate. 

• Allows for sustained wide (> 5 kpc) orbit BHs. 

Correctly modeling the rich orbital dynamics of a black 
hole within its host galaxy can have important consequences 
for its accretion history, duty cycle and observability that 
was previously neglected in simplified ‘advection’ schemes. 
The dynamically complex and more realistic orbits allowed 
by our method will have crucial implications for the early 
growth of SMBHs, which takes p lace in small, rapidl y 
growing galaxies at high redshift (lAvkutalp et al.l l2014h . 
Additionally, understanding the relative importance of dif¬ 
ferent accretion mechanisms throughout a SMBH’s lifetime 
requires the ability to accurately model its dynamics during 


all phases of galaxy evolution, including merger events. 
Implementation of DF routines such as the one presented 
here will improve the ability of cosmological simulations 
to accurately model SMBH accretion, growth and energy 
deposition in the IGM and, therefore, increase the ability of 
simulations to interpret and predict observational results. 
The implementation of this approach in our future cosmo¬ 
logical volume and zoom simulations represents an exciting 
chance to realistically study the growth and merger rate of 
SMBHs across cosmic time. 
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